%% Solve model


% Save parameters
cd   'model'
save paramfile sig_eps_i sig_eps_a betta siggma varphi eta ...
     theta_1 theta_2 theta_3 theta_4 theta_5 phi_pi phi_y pibar ...
     rho_a rho_i pbar_1 pbar_2 pbar_3 pbar_4 pbar_5 ...
     ybar mctildebar pstarbar_1 pstarbar_2 pstarbar_3 pstarbar_4 pstarbar_5 ...
     sbar_1 sbar_2 sbar_3 sbar_4 sbar_5 sbar ...
     markupbar_1 markupbar_2 markupbar_3 markupbar_4 markupbar_5 ...
     psibar_1 psibar_2 psibar_3 psibar_4 psibar_5 ...
     phibar_1 phibar_2 phibar_3 phibar_4 phibar_5 ...
     abar nbar tfpbar ibar i_annbar pi_annbar posmp


% Solve 3rd order approximation around deterministic steady state
order = 3;
dynare model_02_specific.mod nograph -Dalternative_taylor_rule=0
cd ..


% Compute linearized model dynamics around stochastic steady state
type = 'stochastic_steady_state';
[deflect_sss] = compute_deflected_linear_approximation(M_,options_,oo_,type);
[deflect_sss] = compute_deflected_irf(M_,options_,deflect_sss,var_list_,0); % last argument: do plots


% Position of variables
fun_var_pos


% IRFs: MP shock
irf_mp_sss = deflect_sss.irf(1:M_.endo_nbr,:)';

% level IRFs
T = size(irf_mp_sss,1); 
irf_mp_sss_lev = exp( deflect_sss.y' + irf_mp_sss );


% IRFs: price level 
sss_p_1 = exp(deflect_sss.y(pos.pstar_1));
sss_p_2 = exp(deflect_sss.y(pos.pstar_2));
sss_p_3 = exp(deflect_sss.y(pos.pstar_3));
sss_p_4 = exp(deflect_sss.y(pos.pstar_4));
sss_p_5 = exp(deflect_sss.y(pos.pstar_5));
irf_mp_sss_p_1(1,1) = ( (1-theta_1) * irf_mp_sss_lev(1,pos.pstar_1) ^ (1-eta) ...
                         + theta_1  * irf_mp_sss_lev(1,pos.pi) ^ (eta-1) ... 
                                    * sss_p_1 ^ (1-eta) ) ^ (1/(1-eta));
irf_mp_sss_p_2(1,1) = ( (1-theta_2) * irf_mp_sss_lev(1,pos.pstar_2) ^ (1-eta) ...
                         + theta_2  * irf_mp_sss_lev(1,pos.pi) ^ (eta-1) ... 
                                    * sss_p_2 ^ (1-eta) ) ^ (1/(1-eta));
irf_mp_sss_p_3(1,1) = ( (1-theta_3) * irf_mp_sss_lev(1,pos.pstar_3) ^ (1-eta) ...
                         + theta_3  * irf_mp_sss_lev(1,pos.pi) ^ (eta-1) ... 
                                    * sss_p_3 ^ (1-eta) ) ^ (1/(1-eta));
irf_mp_sss_p_4(1,1) = ( (1-theta_4) * irf_mp_sss_lev(1,pos.pstar_4) ^ (1-eta) ...
                         + theta_4  * irf_mp_sss_lev(1,pos.pi) ^ (eta-1) ... 
                                    * sss_p_4 ^ (1-eta) ) ^ (1/(1-eta));
irf_mp_sss_p_5(1,1) = ( (1-theta_5) * irf_mp_sss_lev(1,pos.pstar_5) ^ (1-eta) ...
                         + theta_5  * irf_mp_sss_lev(1,pos.pi) ^ (eta-1) ... 
                                    * sss_p_5 ^ (1-eta) ) ^ (1/(1-eta));
for t=2:T
    irf_mp_sss_p_1(t,1) = ( (1-theta_1) * irf_mp_sss_lev(t,pos.pstar_1) ^ (1-eta) ...
                             + theta_1  * irf_mp_sss_lev(t,pos.pi) ^ (eta-1) ... 
                                        * irf_mp_sss_p_1(t-1) ^ (1-eta) ) ^ (1/(1-eta));
    irf_mp_sss_p_2(t,1) = ( (1-theta_2) * irf_mp_sss_lev(t,pos.pstar_2) ^ (1-eta) ...
                             + theta_2  * irf_mp_sss_lev(t,pos.pi) ^ (eta-1) ... 
                                        * irf_mp_sss_p_2(t-1) ^ (1-eta) ) ^ (1/(1-eta));
    irf_mp_sss_p_3(t,1) = ( (1-theta_3) * irf_mp_sss_lev(t,pos.pstar_3) ^ (1-eta) ...
                             + theta_3  * irf_mp_sss_lev(t,pos.pi) ^ (eta-1) ... 
                                        * irf_mp_sss_p_3(t-1) ^ (1-eta) ) ^ (1/(1-eta));
    irf_mp_sss_p_4(t,1) = ( (1-theta_4) * irf_mp_sss_lev(t,pos.pstar_4) ^ (1-eta) ...
                             + theta_4  * irf_mp_sss_lev(t,pos.pi) ^ (eta-1) ... 
                                        * irf_mp_sss_p_4(t-1) ^ (1-eta) ) ^ (1/(1-eta));
    irf_mp_sss_p_5(t,1) = ( (1-theta_5) * irf_mp_sss_lev(t,pos.pstar_5) ^ (1-eta) ...
                             + theta_5  * irf_mp_sss_lev(t,pos.pi) ^ (eta-1) ... 
                                        * irf_mp_sss_p_5(t-1) ^ (1-eta) ) ^ (1/(1-eta));
end
irf_mp_sss(:,pos.p_1) = log(irf_mp_sss_p_1) - log(sss_p_1);
irf_mp_sss(:,pos.p_2) = log(irf_mp_sss_p_2) - log(sss_p_2);
irf_mp_sss(:,pos.p_3) = log(irf_mp_sss_p_3) - log(sss_p_3);
irf_mp_sss(:,pos.p_4) = log(irf_mp_sss_p_4) - log(sss_p_4);
irf_mp_sss(:,pos.p_5) = log(irf_mp_sss_p_5) - log(sss_p_5);


% IRFs: markup level
sss_markup_1 = exp(deflect_sss.y(pos.pstar_1)) ./ ( exp(deflect_sss.y(pos.mctilde)) .* exp(deflect_sss.y(pos.pstar_1)) .^ (-eta*varphi) );
sss_markup_2 = exp(deflect_sss.y(pos.pstar_2)) ./ ( exp(deflect_sss.y(pos.mctilde)) .* exp(deflect_sss.y(pos.pstar_2)) .^ (-eta*varphi) );
sss_markup_3 = exp(deflect_sss.y(pos.pstar_3)) ./ ( exp(deflect_sss.y(pos.mctilde)) .* exp(deflect_sss.y(pos.pstar_3)) .^ (-eta*varphi) );
sss_markup_4 = exp(deflect_sss.y(pos.pstar_4)) ./ ( exp(deflect_sss.y(pos.mctilde)) .* exp(deflect_sss.y(pos.pstar_4)) .^ (-eta*varphi) );
sss_markup_5 = exp(deflect_sss.y(pos.pstar_5)) ./ ( exp(deflect_sss.y(pos.mctilde)) .* exp(deflect_sss.y(pos.pstar_5)) .^ (-eta*varphi) );
irf_mp_sss_markup_1 = irf_mp_sss_p_1 ./ ( irf_mp_sss_lev(:,pos.mctilde) .* irf_mp_sss_p_1.^(-eta*varphi) );
irf_mp_sss_markup_2 = irf_mp_sss_p_2 ./ ( irf_mp_sss_lev(:,pos.mctilde) .* irf_mp_sss_p_2.^(-eta*varphi) );
irf_mp_sss_markup_3 = irf_mp_sss_p_3 ./ ( irf_mp_sss_lev(:,pos.mctilde) .* irf_mp_sss_p_3.^(-eta*varphi) );
irf_mp_sss_markup_4 = irf_mp_sss_p_4 ./ ( irf_mp_sss_lev(:,pos.mctilde) .* irf_mp_sss_p_4.^(-eta*varphi) );
irf_mp_sss_markup_5 = irf_mp_sss_p_5 ./ ( irf_mp_sss_lev(:,pos.mctilde) .* irf_mp_sss_p_5.^(-eta*varphi) );
irf_mp_sss(:,pos.markup_1) = log(irf_mp_sss_markup_1) - log(sss_markup_1);
irf_mp_sss(:,pos.markup_2) = log(irf_mp_sss_markup_2) - log(sss_markup_2);
irf_mp_sss(:,pos.markup_3) = log(irf_mp_sss_markup_3) - log(sss_markup_3);
irf_mp_sss(:,pos.markup_4) = log(irf_mp_sss_markup_4) - log(sss_markup_4);
irf_mp_sss(:,pos.markup_5) = log(irf_mp_sss_markup_5) - log(sss_markup_5);


% IRFs: markup dispersion 
sss_avg_markup = 1/5*(sss_markup_1 + sss_markup_2 + sss_markup_3 + sss_markup_4 + sss_markup_5 );
sss_markupdisp = 1/5*(sss_markup_1 - sss_avg_markup).^2 + 1/5*(sss_markup_2 - sss_avg_markup).^2 + 1/5*(sss_markup_3 - sss_avg_markup).^2 + 1/5*(sss_markup_4 - sss_avg_markup).^2 + 1/5*(sss_markup_5 - sss_avg_markup).^2;
irf_mp_sss_avg_markup = 1/5*(irf_mp_sss_markup_1 + irf_mp_sss_markup_2 + irf_mp_sss_markup_3 + irf_mp_sss_markup_4 + irf_mp_sss_markup_5 );
irf_mp_sss_markupdisp = 1/5*(irf_mp_sss_markup_1 - irf_mp_sss_avg_markup).^2 + 1/5*(irf_mp_sss_markup_2 - irf_mp_sss_avg_markup).^2 + 1/5*(irf_mp_sss_markup_3 - irf_mp_sss_avg_markup).^2 + 1/5*(irf_mp_sss_markup_4 - irf_mp_sss_avg_markup).^2 + 1/5*(irf_mp_sss_markup_5 - irf_mp_sss_avg_markup).^2;
pos.markupdisp = M_.endo_nbr+1;
irf_mp_sss(:,pos.markupdisp) = irf_mp_sss_markupdisp - sss_markupdisp;
irf_mp_sss(:,pos.tfp) = - eta / 2 * irf_mp_sss(:,pos.markupdisp);


% IRFs: labor 
irf_mp_sss(:,pos.n) = irf_mp_sss(:,pos.y) + eta/2 * irf_mp_sss(:,pos.markupdisp);


% IRFs: aggregate markup
pos.agg_markup = M_.endo_nbr+2;
irf_mp_sss(:,pos.agg_markup) = -irf_mp_sss(:,pos.mctilde)+irf_mp_sss(:,pos.tfp);


